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A  RECURSIVE  ESTIMATION  ALGORITHM 
FOR  IDENTIFICATION  OF  DYNAMIC  SYSTEMS 

ABSTRACT 

This  report  deals  with  the  identification  of  unknown  parameters 
in  dynamic  systems.  In  modeling  a  physical  system,  the  problem  of 
identifying  the  dynamics  of  a  system  are  often  encountered.  The 
developed  algorithm  provides  a  tool  to  model  all  or  parts  of  a  dynamic 
system  using  input-output  data  sets  from  a  real  system.  The  methodol¬ 
ogy  and  techniques  of  this  algorithm  arc  based  upon  linear  recursive 
estimation  theory.  The  theoretical  foundation  and  the  pragmatics  of 
utilizing  the  ensemble  data  to  estimate  the  unknown  parameters  are 
discussed  at  length  in  the  development  of  the  algorithm. 

As  an  application  of  the  algorithm,  experimental  data  from  a 
man-in-the-loop  simulation  is  used  to  estimate  the  parameters  of  a 
single  axis  model  of  the  gunner.  The  tracking  response  of  the  gunner 
model  compare  favorably  with  data  obtained  from  the  simulation.  The 
differences  in  tracking  responses  are  attributed  to  not  including 
human  randomness  in  the  model. 
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A  RECURSIVE  ESTIMATION  ALGORITHM 
FOR  IDENTIFICATION  OF  DYNAMIC  SYSTEMS 


1 .  INTRODUCTION 


This  report  is  a  discussion  and  documentation  of  an  algorithm  for 
identification  of  dynamic  systems.  The  proposed  algorithm  is  closely 
related  to  other  sequential  estimation  algorithms  proposed  by  Kalma^^ 
Mayne f^Deutch^  and  others  over  the  past  decade.  The  most  signifi¬ 
cant  and  unique  feature  of  the  algorithm  described  here  is  the  extent 
to  which  computational  efficiency  and  accuracy  are  achieved  and  assur¬ 
ed.  Particular  attention  is  given  to  handling  of  ill-conditioned 
computations,  and  non-stationary  models. 


The  algorithm  described  here  was  first  implemented  on  a  small 
computer  for  a  real  time  application  and  subsequently  has  continued 
to  be  refined  for  small,  highly  interactive  computers.  The  interactive 
features  of  the  algorithm  have  proven  to  be  invaluable  in  all  stages 
of  an  estimation  project.  The  analyst  may  interact  early  in  an  applica¬ 
tion  to  tune  the  free  parameters  of  the  algorithm  and  throughout  pro¬ 
duction  of  estimates  to  control  computing  time,  convergence  of  the 
estimates,  and  computer  output. 


The  viewpoint  taken  in  developing  and  listing  the  algorithm  des¬ 
cribed  here  was  that  the  data  upon  which  the  estimates  are  based 
arrive  one  at  a  time  in  a  discrete  unedited  form.  The  capacity  to 
store  old  data  and  the  computer  time  are  both  very  limited  resources. 
The  objective  of  the  algorithm  is  to  compute  a  new  model  from  each  new 
data  set  utilizing  a  compacted  form  of  the  data  history.  The  extent 
to  which  the  history  is  utilized  is  controlled  by  the  algorithm  in 
such  a  way  as  to  minimize  the  prediction  error  for  the  model. 


The  remainder  of  this  report  is  a  discussion  of  the  theoretical 
foundations  of  the  algorithm  and  the  pragmatics  of  data  collection  and 
model  validation. 
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2. 


THEORY 


2.1  Statistical  Model. 

The  basic  model  of  the  random  process  to  be  observed  is 

Y  =X  M  V  (NX  1)  2.1.1 

Where  Y  is  the  vector  of  observations,  X  (N  X  P)  is  the  matrix  of 
controlled  variables,  6  (P  X  1)  is  the  vector  of  parameters  to  be 
estimated,  and  V  (N  X  1)  is  the  vector  of  observation  errors.  The 
observation  errors  are  presumed  to  have  well  known  first  and  second 
moments  defined  by 

E (V)  =  0  2.1.2 

E(VVT)  =i|j  2.1.3 

A 

We  wish  to  compute  an  estimate,  8_,  of  such  that 

A 

B_  =  K  Y  2.1.4 

A 

and  E (B)  =  i  (unbiased) 

This  implies  directly  that 

XX  -  I_  =  0  for  aii  B.  2.1.5 

The  estimation  problem  may  be  formulated  as  a  minimum  variance  problem 
where  J(  must  be  found  such  that 

Var (£)  =  Tr  E  [ ( B - B) ( B -  8) T]  2.1.6 

in  minimized,  subject  to  Equation  2.1.5.  The  minimization  of  Equation 
2.1.6  leads  directly  to 

T  -1  -IT  -1 

K  =  (X  \p  X)  V£  2.1.7 

2 . 2  Model  of  the  Observation  Errors. 

In  the  previous  section  it  was  presumed  that  the  covariance 
of  the  observation  errors  was  known.  We  will  now  develop  a  statistical 
model  of  the  observation  errors  which  produces  a  sequence  of  covariance 
matrices,  i^,  where  i  is  conceptually  an  index  on  time.  The  model  will 
represent  an  acknowledgement  that  observation  errors  increase  with  the 
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age  of  the  observation. 


t 


i+1 


O 

0 


4-aJ^ 


2.2.1 


In  this  model  variances  occur  at  the  value,  —  and  increase  by  a  factor 

1  .  “ 

— —  in  each  time  increment.  The  variances  are  thus  increasing  exponent¬ 
ially  with  age  such  that  diagonal  elements,  t  units  of  time  in  age, 

,  ,,  „  -t  o ^ 

are  equal  to  (1-a)  —  . 

It  is  convenient  in  the  development  which  follows  to  repre¬ 
sent  the  error  variance  model  as 

.2 


ld.i  ■  0  2i*i 


2.2.2 


where 


S-i+1 


I  i 

a  | 


T--r:o: 

!  W*1 


2.2.3 


Q„  -  1/a 


Using  Equation  2.2.2,  the  estimation  gain  matrix  K,  reduces  to 

K  =  (X.  .oT1.  X.  .)_1X.  .  QT1. 

—  — l+l-^-i+l  —l+l  —l+l  -^i+I 


2.2.4 


Where  X^+,  same  as  X.  of  Equation  2.1.7.  The  subscript  of  X^+j 

affords  a  means  of  writing  a  recursive  definition 


X.  , 
—l+l 


L+l 


X. 

—1 


2.2.5 


where  x.  .  is  the  vector  of  most  recent  values  of  the  controlled  variables. 
Likewise  the  observations  may  be  defined  recursively  as 


Ii+1 


fi+1 

Y. 

—l 


2.2.6 


Using  Equations  2.2.4,  2.2.5,  and  2.2.6  the  estimation  equation 
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2.2.7 


(Equation  2.1.4)  now  becomes 

i.i  ■  eli  fill  W1  eli  ail  ii.i 

A 

Further  use  of  the  stepwise  recursion  for  8  leads  to 

(ax.  .  xT  .  +  (l-a)xTo7*X. )A. 8  =  ax.  .  (y.  ,-xT  ,80  2.2.8 

v  —l+l  —l+l  v  '— 1-%  — i  i—  — l+l i+I  — i+l— i' 


where 


-  £in  -  Si 


Equation  2.2.8  represents  the  starting  point  for  the  develop¬ 
ment  of  the  recursive  estimation  algorithm.  The  assumption  involved 
in  the  use  of  this  equation  are 

1)  Observation  model  (Equation  2.1.1) 

2)  Linear  Estimation  (Equation  2.1.4) 

3)  Unbiased  Estimation  (Equation  2.1.5) 

4)  Minimum  Variance  Estimation  (Equation  2.1.6) 

5)  Known  Error  Moments  (Equation  2.1.2  and  2.1.3) 

6)  Exponential  Error  Variances  (Equation  2.2.1) 

Of  particular  note  here  is  the  fact  that  no  assumptions  are  made  re¬ 
garding  distribution  functions.  Such  assumptions  may  be  delayed  until 

/s 

the  estimates,  and  their  variances  must  be  interpreted  statistically. 

2.3  Statistical  Interpretation. 

It  was  shown  in  the  previous  development  that  the  recursive 
estimation  scheme  (Equation  2.2.8)  produces  an  unbiased,  minimum  var- 

A  A 

iance  estimate  of  the  parameters,  8.  The  estimator,  £,  is  itself  a 
random  variable  thus  a  statistical  measure  of  its  closeness  to  the  true 

A 

parameters,  3,  is  needed.  Such  a  measure  is  the  covariance  matrix. 


P. 


P.  =  E 

—l+l 


['Li  -  Si.pLi 


-  «wT] 


2.3.1 


Using  Equation  2.2.7  and  Equation  2.1.1  the  covariance  matrix  may  be 
written  as 

P.  .  =  o2(xT  .  QT*  X.  .r1  2.3.2 

—l+l  '■—l+l  -*1+1  —l+l 
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2.3.3 


Once  again  a  recursive  formulation  may  be  employed 
Li+1  =  °2  Cca.+1  xj+1  +  a2(l-a)  P'1)"1 


or  equivalently  by  the  matrix  inversion  lemma 


(1) 


P.  .  =  t— — [P .  -  a  P.x.  .  (axT  ,P.x.  ,  +  (l-a)a2)_1xT  .P.l  2.3.4 

—l+l  l-al— i  —l—i+l  v  —l+l— l—i+l  v  7  7  — i+l— iJ 

The  preceding  result  illustrates  how  the  covariance  matrix  may  be  pro¬ 
pagated  forward  in  time  without  the  inverse  indicated  in  Equation 
2.3.3.  In  fact  only  scalar  inversion  is  required  in  Equation  2.3.4. 

In  the  development  of  Equations  2.3.3  and  2.3.4  it  was  pre¬ 
sumed  that  the  necessary  matrix  inverses  always  exist.  In  particular 
T  -1  -1 

the  existence  of()L+^  was  never  questioned.  It  is  quite 

possible  in  the  progression  of  this  matrix  that  it  becomes  singular 

for  all  practical  purposes.  For  example,  if  the  data  became  constant 

(i.e.  x_.  =  x..  for  all  i),  then  the  recallable  history  will  look 

i+l  l  T  -1 

like  a  single  sample.  The  rank  of  the  matrix,  )L  fhen 

one. 

In  the  case  of  singular  or  very  poorly  conditioned  solutions 
of  Equation  2.2.8,  we  may  wish  to  employ  a  generalized  inverse  which 

A 

would  minimize  some  suitable  norm  of  A. 3*  Thus  with  an  incomplete  data 
set  we  would  hold  the  parameter  estimates,  as  constant  as  possible. 
The  covariance  matrix  in  this  case  may  be  found  by  first  rearranging 
Equation  2.2.7  and  substituting  Equation  2.1.1. 

«I*1  ill !  wii.!  ■  *1.1  Sill  4*1  ®i.l  *  *1*1  2i'*A*l 


Now  employing  a  generalized  inverse  leads  to 

k.l  -  4*1  ■  '*1*1  Siil  *i*l’T*I*l  ^1  Vi.! 


2.3.5 


2.3.6 


where,  T,  indicates  generalized  inverse.  The  covariance  matrix  follows 
from  its  definition 

T  -1  T  -1  T  T  7  T  -1 

_  2  ,VT  ^-1  v  ,T  X  ,Q.  ,X.  ,[(x!  ,Q.  .X.  ,)  1  =o  (x!  ,0.  ,X.  J 

P.  =  a  (X.  .0.  ,  X.  ,)  — i+1-H+l-i+l 1  — 1+I-+1+I— i+l7  J  — l+l^i+l— i+l7 

—i+l  — i+l-H+l— i+l 

2.3.7 
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The  interpretation  of  the  covariance  matrix,  is  closely  connected 

to  the  interpretation  of  the  error  covariance  matrix,  ^ .  If  the 
latter  represents  the  second  moment  of  a  Gaussian  distribution,  then 
may  be  interpreted  as  the  second  moment  of  a  Gaussian  distribu¬ 
tion  of  the  estimate  3.  The  covariance  matrix,  P.  , ,  defines  the 

—  *  —l+l’ 

variance  ellipsoids  analogous  to  the  familiar  a,  2a,  3a,  etc.  limits 
of  a  Gaussian  distribution  in  one  dimension. 

2.4  Computational  Considerations. 

From  a  computational  point  of  view,  the  recursive  estima¬ 
tion  scheme.  Equation  2.2.8  is  a  linear  algebraic  equation  of  the  form 


A.  ,  Z.  =  B.  . 

—i+l  —i+l  —i+l 

2.4.1 

where  A.  =  (a  x.  .x!  .  +  (1-a)  xT  Q.^  X.) 

—i+l  v  —i+l— i+l  v  1  —i  — l  -iJ 

2.4.2 

Z.  .  =  A. 3 

—i+l  i— 

2.4.3 

B.  .  5  a  x.  ,  (y.  ,  -  xT  .  3- ) 

—i+l  — i+lv/i+l  —i+l  — \J 

2.4.4 

The  matrix.  A,  is  symmetric  and  is  of  uncertain  rank.  Both  A  and  IJ  are 
functions  of  the  parameter,  a,  which  hereafter  will  be  viewed  as  a  para¬ 
meter  which  may  be  reset  at  each  time.  Thus  any  information  generated 

A 

in  computing  the  estimator,  may  also  be  used  to  determine  the 

covariance  of  the  observation  errors.  Thus  the  statistics  of  the 
observation  errors  may  be  subjected  to  adaptation.  Since  it  is  intend- 

A 

ed  here  that  3^+^  be  used  as  a  predictor  in  the  interval  (i+1,  i+2),  the 
most  appropriate  value  of  a  is  that  which  minimizes  the  single  step 

A 

prediction  error.  If  such  an  a  can  be  found  for  3-,  it  will  simply  be 
used  to  compute  The  single  step  prediction  error,  E^,  is  defined 


by 


E  *  (y.  .  -  xT  .  3.) 

r  w  i+l  —l+l  — iJ 


2.4.5 


The  minimum  of  the  squared  prediction  error,  a  ,  will  be  approximated 
by  a  one  step  quadient  search  over  the  parameter,  a 
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where  a.,,  -  -  c  ^ 

2.4.6 

2.4.7 

2.4.8 

Note  that  in  this  adjustment  of  a  the  only  matrix  inversion  required 

A 

is  the  same  one  completed  in  the  computation  of  By  such  a  scheme 
a  new  value  of  a  may  be  found  with  very  little  computation  over  and 

A 

above  that  required  to  compute  3- 

A 

In  computing  the  estimator,  it  must  be  recognized  that 

A^^,  may  become  nearly  singular.  In  such  degenerate  cases  the  computa¬ 
tion  of  the  estimator  is  very  poorly  conditioned.  In  order  to  have  a 
measure  of  approach  to  singularity,  a  recursive  bound  for  the  eigen¬ 
values  of  is  computed.  Using  Equation  2.4.2  and  the  Rayleigh 

Quotient  it  can  be  shown  that 


<Va.  ,  i  ‘'-“HVa. 

—l+l  —l 

2.4.9 

(X  )  <  (1-a) (X  ).  +  axT  .  x.  . 

v  n'A.  ,  —  v  n  A.  — l+l  — l+l 

2.4.10 

where  A ^  and  X are  the  smallest  and  largest  eigenvalues  respectively. 
((A^a  should  be  read  ’’the  ith  eigenvalue  of  A”).  By  computing  these 
bounds,  the  condition  of  A^+1  can  be  approximated  and  used  to  trigger 
a  generalized  inversion  whenever  necessary. 

2 . 5  Identification  of  Discrete  Systems 

The  recursive  estimation  technique  developed  in  the  previous 
section  will  now  be  used  as  a  basis  for  identification  of  dynamical 
systems.  As  an  example  suppose  the  system  to  be  identified  is  modelled 
as  a  linear  discrete  system  of  the  form, 

Xi+1  =  £  (i,  i  +  1)  X.+  B  U.  2.5.1 
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where  )L  =  state  at  time  i 

£(i,  i+1)  =  state  transition  matrix  from  i  to  i+1. 

As  an  example  consider  a  single  input,  single  output  system. 


In  such  a  system  the  state  vector  will  consist  of 
u. 


X. 

—l 


E-1  x. 

p-2  1 
E  x. 
1 


2.5.2 


and  the  input  forcing  vector  may  be 


U.  = 

—l 


u. 

—l 

p"1 

E  u. 

l 

p-2 

E  Ui 


where  Enx.  =  x.  =  Z+nx. 

l  l+n  i 


2.5.3 


(En  =  E(E(E( - (•)( - );  n  applications  of  E) 

n  times 

In  this  example  suppose  the  Z- transfer  function  was 


Gm  -  eZ5  +  fz2 

“43  2 

Z  -aZ  -bZ  -cZ-d 


2.5.4 


The  state  variable  representation  of  such  a  system  is 


X.  . 
—l+l 


a 

b 

c 

d  ~ 

e 

f~ 

1 

0 

0 

0 

0 

0 

_ 

X.  + 

0 

1 

0 

0 

—1 

0 

0 

0 

0 

1 

0 

0 

0 

U. 

—l 


2.5.5 
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where 


The  objective  of  an  identification  algorithm  for  such  a 
situation  is  the  estimation  of  the  values  of  a,  b,  c,  d,  e,  and  f. 
The  recursive  estimation  algorithm  described  in  the  previous  section 
is  directly  applicable  after  the  following  transformations  are  made 

Discrete  System 
Variables 


x. 

1 

txi-l'  xi_2’  xi-3’  Xi-4 ’ 

Ui-1*  Ui-2] 

[a,b,c,d,e,f] 

The  viewpoint  employed  here  is  that  the  most  recent  output  observation 
is  subject  to  error,  e.  There  is  an  underlying  assumption  that  the 
other  variables,  X^,  are  error  free.  To  the  extent  that  such  an  assump¬ 
tion  is  not  valid,  these  errors  may  be  thought  of  as  included  in  e 
by  the  following 

c .  =  Error  in  x.  +  6^  , •  [Error  in  X.  ]  2.5.6 

l  i  i-l  — i - 1 

A  more  formal  representation  of  these  errors  is  not  tractable  at  the 
present  time,  either  by  conventional  regression  analysis  or  by  filtering 
theory. 

2 . 6  Identification  of  Continuous  Systems. 

In  contrast  to  the  previous  section  we  will  now  develop  a 
formulation  of  a  continuous  system  which  will  permit  its  identification 
by  the  proposed  algorithm.  Most  of  the  information  suggested  in  the 


Estimation 

Variables 


yi 


T 

X 

—l 


xi  -i 


xi  -2 


Xi  -3 


u. 

1 


Vi 
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literature  require  the  use  of  derivatives  of  state  variables  as 
fl  41 

observations.  ’  This  is  to  be  avoided  if  possible  because  of  the 
excessive  errors  involved  in  observing  or  approximating  these  derivatives. 

the  system  to  be  considered  here  is  of  the  form, 

X  =  A  X  +  B  U  (Nxl)  2.6.1 

Where  A  and  jl  are  the  parameters  (perhaps  time  dependent)  which  must 
be  estimated.  The  analytic  solution  of  this  equation  over  the  time 
interval  is 

X.  .  MX.  +  R  U. 

—l+l  x  — i - i 

where  £  =  exP  (AAt)  2.6.2 

R  =  (£  -  I)  A-1B 

and,  where  for  the  moment,  A  and  are  treated  as  constants.  Now 
Equation  (2.6.2)  is  of  the  same  form  as  Equation  2.5.1,  but  the  state 
variables  are  not  interrelated  as  in  Equations  2.5.2  and  2.5.3.  Here 
the  state  variables  have  been  carried  over  from  the  continuous  form, 
Equation  (2.6.1). 

The  identification  problem  has  now  been  cast  in  the  form  of  a 
problem  of  estimating  £  and  R.  Using  the  algorithm  described  in  the 
section  on  multivariate  observations,  the  estimates  are  obtained  by 


first  transforming 

Equation  (2.6.2) 

Estimation  Equations 

<£_  and  R 

8 

X.  and  U. 

X. 

—i  —l 

—l 

*d+i 

A 

4 

The  estimates,  represent  an  identification  of  the  discrete  form  of 
the  original  system.  The  conversion  back  to  the  continuous  form  is 
most  readily  accomplished  by 
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1)  Noting  that  £  and  A  have  the  same  eigenvalues,  first 
diagonalize  £ 

£  =  S  DjSf 1  =  e^T  =  SAAV 1 

5a  ■  i  L"  4 

A  =  S  D.S 
— - A— 


2.6.3 

2.6.4 

2.6.5 


where  D,  and  are  diagonal izations  of  <f>  and  A  respectively  and  the 
above  computations  may  involve  complex  arithmetic. 


2)  From  A,  <J>,  and  R;  may  be  recovered 
B  =  A  (4)  -  I) -1R 


2.6.6 


2 . 7  Multivariate  Observations. 

The  derivations  of  the  previous  sections  are  readily  extended 
to  include  more  than  one  observation  at  each  time  step.  Beginning  with 
the  model  of  the  random  process 

y.  =  x.  8-  +  V. 

^-l  —i  —l  —i 


E  (Vj)  =  0 


2.7,1 


E (V . V .  )  =  R. 

—l—i  i 


We  wish  to  estimate,  8-,  such  that 


— l J 


J=0  J  J 


2.7.2 


E(e.)  =  b. 


var  (8.)  is  minimized. 

— l 


The  observation  errors  are  readily  modeled  as 


R.  =  -  r 

—l  a  — 


2.7.3 
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R.  . 
— l-t 


U-ct  )ta 


2.7.4 


The  estimates  are  obtained  in  a  form  similar  to  the  Kalman  Filter  for¬ 
mulation 


—l 


or  in  the  special  case  of 


r  =  a  I, 


it  follows  that 


where 


_  1. 

l  -  a‘ 

.-1  . . 


T  -1  ^ 

X.  R.  Cy1  -  X.  Bj.p 

2.7.5 

Lii"1  xj* 

2.7.6 

2.7.7 

Li  *  “Si  il 

lI  ^ 

2.7.8 

P. 

— l 

s  in  the  range 

<  (1-a)  (X  )n-i  +  a(X  )  T 
—  v  1  v  irQ.  ,  v  n'x.  x. 

-*1+1  —l  — i 

2.7.9 

where  should  be  read  "the  jth^  eigenvalue  of  Q"  and 

more  conservatively 

(1-cr)  CXX)  -1  <  (X  )  -1  <  (l-a)(A  )  -1  +  a  Tr(x.Tx.) 

-*i-l  J  iii  -S-i-1 


2.7.10 


2.8  Single  Input  -  Single  Output  Models. 

Much  of  the  analysis  of  the  previous  section  is  simplified, 
and  the  computations  improved  if  one  is  dealing  with  a  single  input- 
single  output  model. 

Consider  for  instance 

x  +  ax  +  bx  =  cU(t)  2.8.1 
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The  derivatives  may  be  approximated  as  follows 


•  I 


X 


i 


-  +  o  (At2) 

At 


=  (Ax^  +  Axi_1)/(2  At)  +  0  (At  ) 


When  these  approximations  are  used  in  the  original  equation  2.8.1,  the 
result  is  a  finite  difference  equation. 


2  2 
a. A  x.  +  a., (Ax.  +  Ax.  .)  +  a7x.  =  cU.  +  0(At  1 
1  l  2*-  i  i-l'  3  i  i  v 


2.8.2 


or 


(oij+o^Ax^  +  (a0-a^)Ax^  ^  +  a^x^  =  clh  +  O(At^) 


2.8.3 


Noting  that 


,  .  2+aAt  ,.-2 

tal  *  a2>  '  2  Lt 


a3  5  b 


Equation  2.8.3  may  be  converted  to 

a3 


Ax.  =  - 

l 


(ara2)  A  c 

-  X .  +  7 - f  Ax.  +  - — 

a^ot^  i  (a^+a^)  i-i  ai+a2  1 


U.  +  0  (At  )  2.8.4 


where  arty  approximation  in  IT  must  be  good  to  within  0(At  ).  The 
identification  problem  is  now  one  of  estimating  the  coefficients  of 
2.8.4  where  the  following  conversion  is  made 


Discrete  System 
Variables 

Ax. 

l 


[x. ,  Ax.  ,  U. ] 
L  i*  i-i,  iJ 


a. 


Va2 


Va2  ai+a2 


al+a2 


Estimation 

Variables 


yi 


T 

X. 

—1 
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3. 


DATA  HANDLING  ASPECTS  OF  ESTIMATION 
3.1  Discussion. 


In  any  actual  estimation  task  there  are  always  a  variety 
of  problems  which  arise  because  of  the  form  of  the  data.  Where  one 
may  wish  to  have  a  continuous  record  of  the  input-output  behavior 
of  the  system,  in  fact  it  may  be  a  finite  record,  discretized  without 
regard  to  the  requirements  of  an  estimation  algorithm.  Many  of  these 
problems  were  actually  encountered  in  the  application  discussed  in  the 
next  section.  There  the  data  was  available  in  finite  records  where 
each  record  represented  one  pass  of  a  target  for  a  man-machine  system. 
An  example  of  the  input  and  output  data  as  a  function  of  time  is  shown 
below. 


There  were  a  great  many  of  the  data  records  as  shown  above.  No  one 
record  contained  sufficient  information  for  estimation  purposes. 

In  order  to  obtain  a  smooth,  extended  data  record,  each 
individual  record  was  repeated  in  a  reverse  direction,  as  illustrated 
below. 

INPUT 

0  TIME  0 

OUTPUT 

0  TIME  0 
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The  change  to  a  new  data  record  was  always  done  at  time  *  0.  Thus  the 
multiple  data  records  were  concatenated  to  form  one  large  record  without 
discontinuities.  With  data  in  the  form  illustrated  above,  the  recursive 
estimation  algorithm  could  be  expected  to  provide  an  updated  model  for 
each  sample  in  the  data. 

3 . 2  Converging  to  a  Universal  Model. 

The  model  obtained  from  a  recursive  estimation  algorithm  will 
change  from  step-to-step  in  the  estimation.  As  new  information  becomes 
available,  the  estimator  will  alter  the  model  to  minimize  the  estimation 
and  prediction  errors.  No  matter  how  accurate  such  a  model  is,  it 
may  be  of  limited  value  if  it  fails  to  condense  the  data  to  a  smaller 
more  fundamental  set. 


At  the  cost  of  somewhat  greater  prediction  errors  the  estimat¬ 
ed  model  may  be  forced  to  converge  to  a  single  "universal  model"  by 
removing  a  from  Equation  2.2.8  which  leads  to 


(x.  .x.  .T  +  XT  Q^X.)  A. 6  =  x.  . 
1  +  1 — 1  +  1  — 1  — 1  — 1  1—  — 1+1 


^l.l  -  il.lD 


3.2.1 


The  variance,  ,  tends  to  zero  thus  forcing  the  convergence  to  the 

estimates . 


3 . 3  Preliminary  Model  Validation. 

A  number  of  controls  are  available  during  estimation  to  pro¬ 
vide  measures  of  performance  of  the  estimated  model.  At  each  step 
of  the  algorithm  discussed  in  the  previous  section,  the  residual  error, 


y.  -  x.  8. 

'  i  —i—i 

and  the  single  step  prediction  error 


yi+i  -  —i+i  Si 


3.3.1 


3.3.2 


are  available.  Likewise,  the  covariance  matrix,  P.,  is  also  readily 


-i' 


available,  but  perhaps  difficult  to  use  without  the  statistics  of  the 
observation  errors  being  specified.  The  correlation  matrix  based  upon 
P.  is  quite  useful  in  determining  the  extent  to  which  the  model  para¬ 
meters  represent  the  same  thing. 
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Once  a  model  is  expected  to  be  valid  for  a  particular  data 
record,  that  model  may  be  used  to  reconstruct  all  the  other  data 
records.  The  reconstructed  record  may  be  compared  with  the  actual 
record  on  the  basis  of  residual  error  and  prediction  error. 

It  is  emphasized  that  the  validation  discussed  here  are 
preliminary.  It  means  little  to  have  a  model  reconstruct  the  data 
whence  it  came.  This  preliminary  validation  must  be  viewed  as  a 
necessary  but  not  sufficient  step  in  the  validation  process.  The 
ultimate  validation  must  be  based  upon  the  model’s  ability  to  "stand 
along"  in  place  of  the  actual  process.  Stability  constraints  may  be 
required  to  insure  that  the  dynamic  model  is  a  valid  replica  of  the 
real  thing. 


4.  APPLICATION  OF  IDENTIFICATION  ALGORITHMS. 

4.1  Discussion. 

The  methodology  developed  in  the  previous  sections  is 
readily  applied  to  the  very  real  problem  of  understanding  the  behavior 
of  man  in  a  tracking  task.  (See  Appendix  B  for  description  of  the 
system  being  controlled  by  man)  As  a  reasonable  cause-and-effect  model 
of  man  we  take 


d+ctjfi+c^S  =  c*3  (0(t-r)  +  a40(t-r)  + 


a50(t-x)) 


where 


6  =  Manual  output,  perturbation  from  nominal 

0  =  Tracking  error 

r  =  Delay  time 


4.1.1 


Using  the  methods  of  Section  2.8,  this  model  is  readily  discretized 
to  become 


iSi-l  ■  bl  Si-l*b24Si-2  *  Vei-l-T*Vei-2-T*b5Ae,i-3-T*0Wt  > 


4.1.2 


where  T  =  T  At 
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The  model  parameters,  b^,  b2>  — —  -  b 5>  may  be  used  to  recover  the 

original  parameters,  ct^,  a2 -  a^.  One  would  expect  both  sets  of 

parameters  to  be  time  varying  and  thus  require  a  robost  estimator  to 
determine  their  values.  Little  will  be  gained,  however,  if  the  estimated 
parameters  are  all  permitted  to  vary  with  time.  The  resulting  model 
will  have  more  time  functions  than  the  original  data  and  as  such  will 
not  represent  a  condensation  of  the  data.  Such  a  condensation  of  the 
data  is  critical  to  establishing  a  fundamental  understanding  of  the 
behavior  of  man. 

Although  it  is  acknowledged  that  man  is  adaptive,  and  thus 
time  varying  in  his  description,  the  search  for  a  universal,  constant 
parameter,  man  is  worth  pursuing.  The  resulting  model,  if  it  is  valid, 
would  represent  a  system  independent,  fundamental  understanding  of 
the  behavior  of  man. 

4 . 2  The  Search  for  a  Universal  Man  Model. 

The  methods  of  Section  2.4  and  2.8  were  adopted  to  the 
search  for  a  universal  man  based  upon  the  model  formulation  of  Equation 
4.1.2.  Performance  data  in  the  form  of  (60)  pairs  were  used  to  estimate 

the  parameters  b^,  bn, - b^.  The  original  estimates  which  were  time 

varying  were  obtained  using  a  fading  memory  estimator.  After  the  pre¬ 
diction  errors  settled  out  to  values  on  the  order  of  0.05  degrees 
manual  rotational  output,  the  fading  memory  was  removed  and  the  para¬ 
meters  allowed  to  seek  universal  values.  The  resulting  model  develop¬ 
ed  a  somewhat  degraded  performance  to  the  level  of  0.1  degrees  of 
rotational  output  error.  The  universal  dimensionless  coefficients 
turned  out  to  be 

bj  =  -0.0103 

b2  =  0.272 

b3  =  0.0368 

t>4  =  8.39  x  10'4  x  17.453 
t>5  =  0.0135 
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with  a  time  delay  of  0.15  sec  associated  with  the  tracking  error  input. 

When  converted  to  the  continuous  form,  Equation  4.1.1,  the 
coefficients  become 


al 

= 

22.9 

sec  ^ 

a2 

= 

6.49 

-2 

sec 

-2 

a3 

= 

404. 

sec 

a4 

= 

.0672 

sec 

.763 

x  10" 

4.3  Validation  of  the  Universal  Man. 

As  indicated  in  the  previous  section,  the  universal  man 
was  able  to  reproduce  the  data  whence  it  was  created  to  within  0.1 
degrees  of  rotational  output  position.  It  soon  became  apparent  that 
such  a  validation  was  not  good  enough.  The  stability  constraints  were 
not  included  in  the  estimation  algorithm,  and  as  a  result  the  continu¬ 
ous  model  of  man  was  incapable  of  handling  the  control  task  in  a 
complete  system  simulation. 

Further  analysis  of  the  data  used  to  force  man  (i.e.  0(t)) 
indicated  that  the  frequency  components  of  0(t)  were  all  of  lower 
frequency  than  natural  frequency  of  man.^-*  This  suggested  that  the 
n euro -muscular  terms  of  the  estimated  man  were  spurious  (i.e.  b^  and 
b0  were  probably  not  estimated  correctly). 

At  this  point  the  estimated  man  model  becomes  a  point  of 
departure  for  an  orderly  search  for  the  correct  n euro -muscular  response 
of  man. 

Previous  studies  indicate  that  the  eigenvalues  of  man  are 
limited  to  at  most  -10  sec-*.^  Taking  one  eigenvalue  at  -10  sec 
a  second  eigenvalue  on  the  order  of  -580  sec  gave  the  most  reason¬ 
able  overall  response.  On  this  basis  the  continuous  representation 
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of  man  becomes 


otj  =  +590  sec 

-2 

a  2  =  5800.  sec 

_2 

a3  =  4040.  sec 

_2 

a.  =  .0672  sec 

4 

=  .763  x  10'3  sec 

The  model  is  limited  by  the  fact  that  it  is  a  localized  model  which  in 
no  way  models  man's  behavior  in  pursuit.  The  response  to  track  error 
is  a  compensatory  response  or  nulling  response.  In  more  mathematical 
terms  it  is  acknowledged  that  the  man  model  obtained  is  a  perturbation 
model  representing  the  motion  of  man  with  respect  to  some  nominal 
motion.  This  nominal  motion  of  man  was  important  in  the  validation 
of  the  man  model  and  was  approximated  by 

^Nominal  K  N1  4.3.2 

where  =  Angular  rate  of  the  target. 

K  =  Constant  on  the  order  of  1.0. 

4.4  Summary  of  the  Validation  Results. 

The  man  model  represented  by  equations  4.3.1  and  4.3.2  was 
used  in  a  number  of  different  tracking  situations.  These  results  are 
plotted  and  shown  in  the  Appendix  C.  The  purpose  of  this  section  is 
to  provide  a  comparison  of  man's  actual  performance  and  the  predicted 
performance  on  the  basis  of  extremes.  One  of  the  most  important 
performance  criteria  to  be  applied  in  judging  the  performance  of  a 
model  is  the  ability  to  limit  the  extreme  values  of  track  error.  The 
performance  of  the  overall  system  is  closely  related  to  man's  ability 
to  limit  the  extremes  of  track  error.  Three  distinctly  different 
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phases  of  performance  are  available  for  comparison  purposes: 

(1)  Acquisition  phase 

(2)  Firing  phase 

(3)  Crossover  phase 

Table  4.4.1  shows  the  performance  of  the  universal  man  model 
(Equations  4.3.1  and  4.3.2  with  K  =  1.1)  and  the  actual  man  during 
these  phases  of  tracking. 
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TABLE  4.4.1 

COMPARISON  OF  MAN  MODEL  WITH  ACTUAL  MAN** 


MAXIMUM  TRACKING  ERRORS,  MILLIRADIANS 


TARGET 

TARGET  ACQUISITION  PHASE 

FIRING 

PHASE 

CROSSOVER 

PHASE 

IDENTIFICATION 

PARAMETERS  MODEL  MAN 

MODEL 

MAN 

MODEL 

MAN 

MMP1 

A=400m  6.2  10.5 

V=231 . 6  m/sec 

X=450m 

3.5 

3.0 

9.4 

9 

MMP3G2 

A=800m  6.7  13.6 

V=231 . 6  m/sec 

X=450m 

4.9 

5.3 

9.4 

13 

MMP5G2 

A=200m  3.4  5 

V=128 . 6  m/sec 

X=450m 

2.2 

2.5 

3.7 

3.6 

MMP6G1 

A=400m  2.7  10.9 

V=180  m/sec 

X=250m 

3.8 

6.7 

24.2 

20 

★ 

A  E  Altitude,  V  E  Velocity,  X  E  Ground  Distance  at  crossover 

** 

Man  Model  is  Equation  4.3.1  and  4.3.2  with  K  =  1.1 
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4.5  Approximations  to  the  Man  Model . 

In  the  final  man  model,  Equations  4.3.1  and  4.3.2,  some  of 
the  parameters  are  of  little  significance  in  describing  the  behavior 
of  man.  In  particular,  the  smallest  eigenvalue,  -580  sec  ,  is  of 
little  significance  since  its  effect  can  only  be  observed  at  very  high 
frequencies.  Likewise,  the  0  term  has  been  found  experimentally  to 
have  very  little  effect  on  the  performance.  Thus  we  might  consider 
the  following  model: 

6  +  0^6  =  cx30(t-t)  +  a40(x-x)  4.5.1 

5xt  .  .  =  K  A. 

Nominal  1 


where 

al  = 

10  sec  * 

a3  = 

6.97  sec” 

a4 

.0672  sec 

K 

1.0  sec 

A  system  block  diagram  of  this  model  is  shown  in  Figure  4.5.1. 


FIGURE  4.5.1  APPROXIMATED  MAN  MODEL 
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APPENDIX  A 


DOCUMENTATION  OF  ESTIMATION  ALGORITHMS 
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MAIN 


PURPOSE:  Overall  Control  of  the  Estimation  Algorithm 

SPECIFIC  FUNCTIONS:  1)  Call  Ready 

2)  Set-up  Estimation  Cycle 

3)  Track  Eigenvalues 

4)  Up-date  Covariance  Matrix,  D=DEMAT 

5)  If  SSW(7)  then  Full  Inverse  with  no 
Weighting 

6)  If  SSW  (1)  then  Stop. 
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FLAT 


PURPOSE:  Controls 
SPECIFIC  FUNCTIONS 


INPUTS  REQUIRED: 


subroutine  flat 


Complete  Deflation  Unconditionally 

1)  Call  on  Orthogonal ization  Routine 

2)  Set  No.  of  Good  Rows  to  1 

3)  Calls  UDATE  for  new  B 

1)  DEMAT 

2)  x,  y 

3)  BI^  (i.e.  previous  model) 


CALL  ORTHO 
ICT  =  1 

CALL  UDATE 
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SUBROUTINE  EIG 


PURPOSE: 

To  Approximate  the  Eigenvalues  of  the  Inverse  of  the 

Covariance  Matrix. 

SPECIFIC 

FUNCTIONS:  1)  Decompose  DEMAT 

2)  Inverse  Power  Iteration  for  RL1 

3)  Power  Iteration  for  RL2 

4)  r / | SZ | :  Largest  absolute  value  of  SZ_ 

INPUTS : 

1)  DEMAT 

OUTPUTS : 

RL1  -  Smallest  Eigenvalue  of  DEMAT 

RL2  -  Largest  Eigenvalue  of  DEMAT 
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SUBROUTINE  EIG 
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SUBROUTINE  DERNK 


PURPOSE:  Controls  Partial  Deflation  of  DEMAT 

SPECIFIC  FUNCTIONS:  1)  Calls  Orthogonalization 

2)  Isolates  Good  Rows  In  Upper  Part  of 

3)  Call  UDATE  for  New  B 

INPUTS:  1)  DEMAT 

2)  x,  y 

3)  BI  (i.e.  previous  model) 


Z 


SUBROUTINE  DERNK 
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SUBROUTINE  ORTHO 


PURPOSE:  To  multiply  DEMAT  by  a  matrix  which  has  mutually 
orthogonal  rows.  The  first  row  is  x. 

SPECIFIC  FUNCTIONS:  1)  Gram  Schmidt  Orthogonalization  with 

respect  to  x^ 

2)  Multiply  result  of  1)  by  DEMAT. 


SUBROUTINE  ORTHO 


39 


SUBROUTINE  FMI 


PURPOSE:  Controls  Full  Inversion  of  DEMAT 

SPECIFIC  FUNCTIONS:  1)  Decompose  DEMAT 

2)  Compute  Right  Hand  Side 

3)  Call  SUBS 

4)  Compute  New  Model,  B^ 

INPUTS:  1)  DEMAT 

2)  x,  y 

3)  BI 


SUBROUTINE  FMI 
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SUBROUTINE  DECOM 


PURPOSE:  Decompose  a  Symmetric  Positive  Definite  Matrix  Stored 
in  Z_ 

SPECIFIC  STEPS:  Choleski  Decomposition 

INPUTS:  1)  Z 

2)  Rank  of  Z  in  ICT 

OUTPUT:  1)  L  of  Z  =  L  J  stored  in  upper  right 

triangular  part  of  Z_ 


SUBROUTINE  DECOM:  ENTER  WITH  ICT  =  NO.  OF  ROWS  OF 

l  (i.e.  ICT  =  RANK  OF  Z) 
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SUBROUTINE  SUBS 


PURPOSE:  Forward  and  Back  Substitution  in  Solution  of 

★ 

L_  Y  =  B_  By  Forward  Substitution 
T 

=  Y  By  Back  Substitution 

SPECIFIC  FUNCTIONS:  1)  Forward  Substitution 

2)  Back  Substitution 

INPUTS:  1)  in  upper  right  part  of  Z_ 

2)  Right  hand  side  in  SZ 

OUTPUTS:  1)  Solution  in  SZ 

SUBROUTINE  SUBS:  J,T  IS  IN  UPPER  TRIANGLE  OF 

RANK  OF  ICT.  SZ  BEGINS  AS 
HAND  SIDE. 


1  WITH 
RIGHT 
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SUBROUTINE  UDATE 


PURPOSE:  Complete  the  Generalized  Inversion  of  the  Matrix,  Z_, 
Having  Known  Rank,  ICT. 

SPECIFIC  FUNCTIONS:  1)  Save  Z  in  ZSAV 

2)  Form  Symmetric  Positive  Definite  Matrix, 

T 

z  +  zz 

* 

3)  Solve  Z_  X_  =  B  for  by  decomposing  Z_, 

★ 

forming  the  right  side,  and  using  SUBS. 

4)  Update  the  Model,  j} 

5)  j}  BI_  +  ZSAVV1  B*  (Pseudo  Inverse) 

INPUTS:  1)  Z,  ICT  (Matrix  and  Rank) 

2)  x,  y 

3)  BI 

NOTE:  Any  estimation  cycle  which  uses  UDATE  must  have  already  used 
ORTHO. 
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SUBROUTINE  UDATE 
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SUBROUTINE  AIM 


PURPOSE: 

SPECIFIC 


INPUT 


To  Compute  a  New  Weighting  Parameter,  a.  Precess  the 
Model . 

9B 

FUNCTIONS:  1)  Set  up  right  hand  side  of  A  — -  B  . 

Conditioned  upon  whether  A  is  full  rank 
or  not.  9b^ 

2)  Solve  for  by  inverting  A  using  new 

right  hand  side. 

3)  Compute  prediction  error  and  smoothed 
prediction  error. 

4)  Compute  new  a  using  steepest  descent. 

5)  BI^  ■*-  jJ;  Precess  the  model  index. 

1)  x>  y 

2)  XOLD,  YOLD 

3)  Decomposed  DEMAT  and  ZSAV 

4)  previous  model. 

5)  EPAIM,  step  size  in  steepest  descent. 
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SUBROUTINE  AIM 
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SUBROUTINE  OBS 


PURPOSE:  Observe  new  data  and  precess  last  observation 

SPECIFIC  FUNCTIONS:  1)  Save  last  observation 

2)  Save  estimation  coefficients  (forward 
pass  thru  data) . 

3)  Observe  new  data 

4)  Compute  squared  prediction  error  and 
squared  current  error  for  each  tracking 
record . 

5)  Validation  of  estimation  coefficients 
with  other  tracking  records  for  any 
target  type. 
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SUBROUTINE  OBS 


SAVE  LAST  OBSERVATIONS 
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SUBROUTINE  READY 


PURPOSE:  General  Initialization 

SPECIFIC  FUNCTIONS:  1)  Clear  arrays  and  set  indices. 

2)  Position  input  tape. 

3)  Name  and  position  optional  output  tape. 

4)  Input  model  configuration 

5)  Initialize  co-variance  matrix  diagonal 
values 


49 


SUBROUTINE  READY 
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APPENDIX  B 


SYSTEM  DESCRIPTION 
AND 

ANALYTICAL  MODEL 
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B.l  INTRODUCTION 


The  Vulcan  Air  Defense  System  (VADS)  is  an  antiaircraft  artillery 
(AAA)  weapon  that  provides  close-in  defense  against  low-flying  aircraft 
and  can  be  used  against  ground  vehicles  such  as  trucks,  personnel 
carriers,  and  other  lightly  armored  vehicles.  The  VADS  utilizes  a  20- 
mm6  barrel  Gatling  type  gun  having  a  fire  rate  of  3000  rounds  per  minute 
with  preselected  burst  lengths  or  1000  rounds  per  minute  with  burst 
lengths  under  operator's  trigger  control.  The  VADS  is  a  manually  operat¬ 
ed  system  equipped  with  a  range-range  rate  only  radar.  The  gunner  is 
required  to  identify,  acquire  and  track  the  target  and  operate  the 
trigger  for  firing. 

The  fire  control  system  uses  a  two  axes  lead  computing  sight, 
mounted  on  the  gun  trunnion,  to  generate  an  angular  displacement  between 
the  gun  tube  and  the  gunner's  line-of-sight .  The  gunner  must  acquire 
the  target  in  the  sight  optics  and  maintain  the  target  in  the  optical 
sight  reticles  (15jri  and  60tf  concentric  circles)  by  controlling  the 
rate  of  gun  motion  in  azimuth  and  elevation.  The  gun  rates  in  conjunc¬ 
tion  with  the  sight  sensitivity  parameter,  T^,  which  is  approximately 
equal  to  ballistic  time  of  flight  of  the  round,  are  used  by  the  lead 
computing  sight  to  generate  the  lead  angle  between  the  gunners  line-of- 
sight  and  gun  position.  The  fire  control  prediction  assumes  the  target 
is  flying  a  straight  line  course  at  constant  velocity. 

The  radar  provides  only  the  present  range  and  range  rate.  The 
instantaneous  range  and  range-rate  data  are  major  inputs  to  the  current 
generator  which  provides  T^  to  the  lead  computing  sight.  The  radar's 
position  relative  to  the  target  depends  upon  the  gunners  ability  to 
track  the  target.  It  is  positioned  in  angle  by  the  gun  mount  position 
in  azimuth  and  by  the  radar's  own  servo  loops,  which  receive  input 
signals  of  gun  elevation  and  two  components  of  lead  angle,  elevation 
and  traverse.  The  radar  must  stay  within  +4  degrees  of  the  target 
to  maintain  lock-on. 
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B. 2  ANALYTICAL  MODELS 


The  mathematical  representation  of  VADS  is  a  nonlinear,  time  varying 
model  of  the  Vulcan  fire  control  system.  The  VADS  model  is  comprised 
of  submodels  for  the  sight  display,  lead  computing  sight,  sight  current 
generator,  turret  servo  and  gun  dynamics,  hand  controller  and  the  operator. 
A  model  of  the  radar  was  developed  but  not  included  in  the  analytical 
model.  Instead,  a  target  model  is  used  to  provide  the  flight  path 
parameter  inputs  to  the  analytical  model.  The  block  diagram  of  the 
VADS  mathematical  model  is  shown  in  Figure  B-l. 

B.2.1  Sight  Display. 

The  sight  display  is  actually  an  integral  part  of  the  lead 
computing  sight,  however,  it  is  considered  either  as  separate  component 
or  combined  with  the  man  model  in  modeling  the  VADS. 

The  sight  display  is  the  movable  sight  reticle.  It  consists 
of  two  concentric  circles,  15  and  60tf  in  diameter.  The  total  field 
of  view  of  the  XM61  sight  is  approximately  400  mils  and  the  diameter 
of  the  display  aperature  is  4  inches.  Hence,  the  sight  display  is 
modeled  as  a  constant,  10  inches/radian. 


B.2.2  Hand  Controller. 

The  operator's  command  to  the  system  is  transmitted  through 
his  manipulation  of  a  hand  controller  (handlebars).  The  handlebars  are 
rotated  like  steering  an  automobile  for  gun  azimuth  commands  and  are 
rotated  vertically  for  elevation  commands  to  the  gun.  The  deflection 
of  the  handlebars  causes  a  voltage  to  be  applied  to  the  inputs  of  the 
rate  servo  loop,  thereby  driving  the  turret  and  gun.  The  output  voltage 
of  the  hand  controller  is  a  nonlinear  function  of  handlebar  position  and 
rate.  The  differential  equation  relating  the  servo  amplifier  voltage, 

V,  to  handlebar  deflection,  6^g,  is 


V  =  K0  6HB  +  kl  6hb 


B.2.1 


where  KQ  is  the  nonlinear  function  relating  voltage  output  to  displace- 
u 
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ment  of  handlebars  and  is  the  aided  laying  tachometer  gain.  This 
nonlinearity  will  be  discussed  in  more  detail  in  subsequent  paragraphs. 
The  gain  can  be  linearized  to  a  constant  without  degrading  the  system 
performance.  The  aided  laying  tach's  main  function  is  to  provide  an 
output  which  is  related  to  rate  of  change  of  handlebar  deflection  during 
the  acquisition  phase. 

B.2.3  Turret  Servo  Loop. 

The  turret  servo  loops  produce  angular  gun  rates  in  azimuth 
and  elevation  in  response  to  handlebar  position,  i.e.  the  turret  sub¬ 
system  is  a  tyle  1  servo  system.  The  equations  of  motion  for  the  turret 
servo  loops  are  made  up  from  the  gun  dynamics,  operational  amplifiers, 
torquer  motors,  tachometers  and  associated  electronics.  A  simplified 
block  diagram  of  the  uncoupled  servo  loop  for  either  axis  is  shown  in 
Figure  B.2.  The  gain,  K^,  in  the  tach  feedback  is  a  nonlinear  function 
of  handlebar  deflection. 

Assuming  that  all  the  parameters  are  constants,  a  transfer 
function  relating  gun  position,  G,  to  voltage  input,  V,  can  easily  be 
derived,  i.e. 


G 

V 


N(Ke  * 


kaW 


s[ 


R  J  T^m 
in  T  m 


Vke  +  kaW 


-s2  + 


r_Jt 

m  I 


Kt(Ke  +  kakxkf) 


-  S  +  1] 


B.2.2 


where 

K  =  Servo  amplifier  gain 

K^,  =  Motor  torque  constant 

K_  =  Motor  back  emf  constant 
E 

Kj  =  Variable  gain  dependent  upon  handlebar  deflection 
and  tachometer  feedback  resistance. 

=  Tachometer  gain 
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Let 


N 

R 


m 


m 


Gain  ratio 
Motor  resistance 

Combined  motor  inertia  and  load  inertia  referred  to 
motor  shaft 


Electrical  time  constant  of  motor 
RmJT 

VKEm+  KAKlV 


B.2.3 


and  K„ 


K. 


G  +  w 

then,  Equation  B.2.2  can  be  written  as 

G  _  KG 

V  S(t  xrS2  +  xrS  +  1 
in  (j  (j 


B.2.4 


B.2.5 


The  time  constant  of  the  motor,  t  ,  is  0.002  sec.  Since  t  is  small 

’  m  m 

(i.e.  t^  1),  Equation  B.2.5  can  be  reduced  and  approximated  by 


V  "  S(tgS  -  1)  B’2,6 

This  is  the  uncoupled  model  of  the  turret  servo  loop. 

Equation  B.2.3  and  B.2.4  show  that  both  Tg  and  Kg  depend  on 
Kj  which  is  a  function  of  handlebar  position,  <$hb*  The  variable  gain 
Kj  cannot  readily  be  computed.  In  order  to  reduce  modeling  errors,  Tg 
and  KG  are  assumed  to  be  constant  and  all  nonlinearity  between  and 
turret  servo  output  are  modeled  in  the  hand  controller's  nonlinear  gain 
KQ.  This  nonlinearity  is  approximated  from  measurements  taken  from  the 
VADS  hybrid  simulator. 

In  the  development  of  the  man  model,  the  operator  is  assumed 
to  be  operating  the  system  in  a  tracking  only  mode  (nonfiring) ,  and  the 
dynamic  coupling  between  the  gun  axes  is  assumed  to  be  zero.  Hence, 
Equation  B.2.6  is  used  as  the  model  of  turret  servo  and  gun  dynamics. 
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B.2.4  Sight  Current  Generator. 

For  the  purpose  of  this  modeling  effort,  the  sight  ourrent 
generator  is  defined  as  the  subsystem  which  produces  the  sight  sensitivity, 
T^,  for  the  lead  computing  sight.  Neither  a  discussion  on  principles 
of  operation  nor  any  developed  mathematical  model  will  be  presented  in 
this  document. ^ 

Instead  of  using  a  analytical  model  of  the  sight  current 

generator,  a  computer  model  is  used  to  generate  the  theoretically  correct 

for  a  straight  and  level  flight  path.  From  the  known  velocity  and 

position  of  the  target  and  the  ballistic  data  of  the  ammunition,  an 

interactive  process  is  used  to  calculate  the  future  range  at  which  a 

projectile  would  intercept  the  target.  Using  the  computed  values  from 

the  interactions,  the  theoretical  T.,  can  be  computed  continuously  using 

C2")  N 

the  following  equations. 


df 

_  4.  rr 

°p 

1  \  1 

DP 

°p  " 

°F 

1-VdD^os  aJ 

B.2.7 


t^  -  projectile  time-of-f light 

cr  -  angle  between  target  flight  path  and  future  range 
vector 

Dp  -  present  range 
Dp  -  future  range 
V  -  velocity  of  target 


For  a  known  straight  and  level  flight  path,  this  equation  provides 
values  of  which  are  approximately  correct. 
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B.2.5  Lead  Computing  Sight. 

The  lead  computing  sight  is  the  "heart"  of  the  Vulcan  fire 
control  system.  It  performs  the  lead  computation  required  to  point 
the  gun  barrels  ahead  of  the  present  position  of  the  target.  Since  the 
projectile  have  a  finite  time  of  flight,  the  sight  predicts  the  future 
position  of  the  target,  where  the  gun  tube  must  be  aimed,  in  order  for 
the  projectile  to  hit  the  target. 

Its  principle  of  operation  is  based  on  the  precessional  law 
of  a  gimballed  gyroscope.  The  case  of  the  sight,  which  houses  an  electro¬ 
magnet,  is  mounted  so  that  it  moves  with  the  gun.  An  eddy  current  disk 
is  attached  to  the  rotor  of  the  two  axis  gimballed  gyroscope.  The  electro¬ 
magnet  is  energized  by  a  current  signal  generated  in  the  sight  current 
generator.  When  a  moving  target  is  tracked,  eddy  current  forces  are 
generated  in  the  disk,  creating  a  torque  that  causes  the  spin  axis 
of  the  gyro  to  follow  the  gun  with  a  displacement  (lag)  sufficient  to 
provide  the  torque  needed  to  precess  the  gyro  at  the  gun  rate.  The 
magnitude  of  the  electromagnet  current  determines  the  angular  lag 
between  the  gun  and  the  gyroscope  spin  axis  for  a  given  tracking  rate. 

This  angle  is  a  measure  of  the  lead  angle  being  developed.  Actually 
the  reticle  orientation  relative  to  the  gun  differs  from  the  gyro 
orientation  relative  to  the  gun  by  a  constant  ratio,  a.  This  permits 
instantaneous  response  between  the  gun  motion  and  reticle  motion.  An 
additional  quantity,  referred  to  as  superelevation,  is  generated  in  the 
elevation  axis  of  the  lead  computing  sight.  Superelevation  compensates 
for  ballistic  drop  of  the  projectile. 

A  simplified  model  of  the  lead  computing  sight  is  developed 
by  considering  the  motion  of  the  sight,  reticle  and  gun  in  a  single 
plane.  Assume  the  Vulcan  is  tracking  a  moving  target  in  the  horizontal 
plane  as  deplicted  in  Figure  B.3.  From  this  figure  the  following  angular 
relations  are  obtained: 

AG  ‘  AGyro  =  gyro  lag  B'2'8 
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X 


B.2.9 


*G  ‘ 


R 


A1  -  Ar  =  9e 
where  =  gun  line  of  sight 

AGyro  =  gyr°  line  of  Sight 


R 


=  reticle  line  of  sight 


=  target  line  of  sight 

X  =  lead  angle 

0  =  track  error 


If  the  sight  gyroscope  with  the  eddy  current  disk  mounted  on  its  spin 
axis  is  rotating  in  an  electromagnet  field  that  is  fixed  to  the  gun, 
the  precessional  rate  of  the  gyroscope  is  related  to  the  displacement 
between  the  gun  and  the  gyroscope  by 


Ar  =  PI  sin  (- 
Gyro 


A„  -  A„ 

G  Gyro 


B.2.11 


where  P,  I  and  K  are  related  to  parameters  of  the  sight  current 
15  K  1 

generator  and  the  lead  computing  sight.  If  PI  =  sr  and  Kc  =  (1  +  a), 

N 

the  small  angle  approximation  of  Equation  B.2.11  becomes 


AG  AGyro 


Gyro  (1  +  a)T. 


B.2.12 


N 


The  reticle  is  controlled  by  the  sight  gyroscope  but  their  lines-of- 
sight  relative  to  the  gun  differ  by  a  fractional  part  of  the  lead  angle, 
oX,  i.e. 


R 


A  =  oX 

Gyro 


B.2.13 
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Substituting  Equation  B.2.13  into  Equation  B.2.9,  the  angular  relations 
between  the  gun  and  gyro  is 

AG  -  Vo  ■  C1  *  0)  X  B-2-14 

The  right  side  of  this  equation  is  the  gyro  lag  relative  to  the  gun.  The 
time  derivative  of  this  equation  is 


*G  -  Vo  =  tl  *  o)  X  B.2.1S 

Solving  Equation  B.2.12,  B.2.14  and  B.2.15  in  terms  of  X  and  A  ,  yields 

S 

a  differential  equation  for  the  lead  angle 


B.2.16 


This  equation  is  only  valid  for  computing  the  lead  angle  in  a  single 
plane. 


In  reality,  the  Vulcan  fire  control  system  is  solving  for 
two  components  of  lead  angle,  traverse  and  elevation.  Since  the  gun, 
reticle  and  sight  gyro  are  not  in  the  same  coordinate  system,  trans¬ 
formulations  are  required  to  project  the  respective  vectors  into  the 
same  coordinate  system.  Neglecting  superelevation  and  using  small  angle 
approximations,  the  dynamic  equations  of  the  sight  in  traverse,  and 
elevation,  A_,  can  be  approximated  by 

C 


Ag  cos  Eg 
1+0 


*T 

Cl  +  a) tn 


B.2.17 


•  _  EG  XE 

aE  (1  +  a)  "  (1  +  a)T, 


N 


B.2.18 


where  Eg  is  the  gun  position  in  elevation.  Although,  the  term,  cos  Eg, 
is  a  liberal  approximation  of  the  required  trnasformulations,  the  lead 
angles  obtained  from  Equations  B.2.17  and  B.2.18  compare  favorable  with 
the  results  obtained  from  experiments  using  the  Vulcan  sight  and  more 
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elaborate  analytical  models  that  take  into  account  the  coordinate 
Cl  31 

transformations.  ' 
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Figure  B.2  Block  Diagram  of  Uncoupled  Turrent  Servo  for  Either  Axes. 


- ►•VELOCITY  OF  TARGET 


Figure  B.3  Single  Axis  Lead  Computing  Sight  Geometry. 
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APPENDIX  C 


TIME  HISTORIES  OF  TRACKING  ERROR 
FROM  VADS  SIMULATION  AND  DEVELOPED  MAN  MODEL 

The  figures  in  this  appendix  show  the  tracking  responses  of  the 
actual  man  and  the  developed  man  model.  The  results  are  shown  for 
four  different  target  conditions.  The  target  is  flying  at  a  constant 
velocity  and  altitude  on  a  straight  fly-by  course. 

The  plots  for  the  actual  man  are  the  averages  of  a  trained  gunner 
tracking  the  same  target  ten  times. 

The  data  were  obtained  from  a  man  operated  simulation  of  the 
Vulcan  Air  Defense  System's  (VADS)  fire  control  system.  The  different 
target  fly-bys  were  designed  to  provide  the  gunner  with  tracking 
tasks  ranging  from  easy  to  difficult. 

The  tracking  response  of  the  model  approaches  the  magnitude  but 
not  the  frequency  of  the  actual  gunner  tracking. 

One  of  man's  limitations  that  tends  to  degrade  system  performance 
is  his  randomness.  This  effect  is  not  included  in  the  developed  man 
model.  The  recursive  estimated  algorithm  smoothed  out  this  randomness 
in  estimated  the  parameters  of  the  model. 
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TRACKING  ERROR,  MILLIRADIANS 


VELOCITY:  2316  m/sec 

ALTITUDE:  400  meters 


Figure  C.l  Tracking  Errors  for  Pass  No.  MMP1. 
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TRACKING  ERROR,  MIlllRADIANS 


TARGET  PARAMETERS 


VELOCITY:  231.6  m/sec 

ALTITUDE:  800  meters 

CROSSOVER  RANGE:  450  meters 


Figure  C.2  Tracking  Errors  for  Pass  No.  MMP3. 
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VELOCITY^  128.6  meters/sec 
ALTITUDE :  200.00  meters 

CROSSOVER  RANGE:  450  meters 


SNviavaimw  'aoaaa  9Ni>ovai 
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Figure  C.3  Tracking  Errors  for  Pass  No  MMP5. 


TRACKING  ERROR,  MILLIRADIANS 


TARGET  PARAMETERS 


71 


Next 


Page  is  blank. 


DISTRIBUTION  LIST 


No.  of 
Copies 

Organization 

No.  of 
Copies 

Organization 

12 

Commander 

Defense  Documentation  Center 
ATTN:  TIPCR 

Cameron  Station 

Alexandria,  VA  22314 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCRD-E 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCCP 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN :  AMCRD-G 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCDL 

5001  Eisenhower  Avenue 
Alexandria,  VA  22314 

1 

Commander 

US  Army  Materiel  Command 
ATTN :  AMCRD-M 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCMA 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCRD-R 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCGS-L 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCRD-P 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCPA-S 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCRD-W 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN:  AMCQA 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1 

Commander 

US  Army  Materiel  Command 
ATTN :  AMCRD-T 

5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

73 


DISTRIBUTION  LIST 


No.  of  No.  of 

Copies  Organization  Copies  Organization 


1  Commander 

US  Army  Materiel  Command 
ATTN:  AMCSA-PM-MBT 
5001  Eisenhower  Avenue 
Alexandria,  VA  22333 

1  Commander 

US  Army  Materiel  Command 
ATTN:  AMCPM-SA 
Picatinny  Arsenal 
Dover,  NJ  07801 

1  Commander 

US  Army  Armament  Command 

ATTN:  AMSAR-SA 

Rock  Island,  IL  61201 

1  Commander 

Harry  Diamond  Laboratories 
ATTN:  Sys  Anal  Ofc 
Washington,  D.C.  20438 

1  Commander 

Frankford  Arsenal 
Philadelphia,  PA  19137 

1  Chief,  Analytical  Sciences 
Office 

USA  Biological  Defense 
Research  Laboratory 
ATTN:  AMXBL-AS 
Dugway,  UT  84022 

1  Commander 

US  Army  Aviation  Systems 
Command 

ATTN:  AMSAV-R-X  (Mr.  Hollis) 
PO  Box  209 

St.  Louis,  MO  63166 
1  Commander 

US  Army  Electronics  Command 

ATTN:  AMSEL-PL 

Fort  Monmouth,  NJ  07703 


1  Commander 

US  Army  Safeguard  Evaluation 
Agency 

ATTN :  SSEA-T 

White  Sands  Missile  Range, 

NM  88002 

1  Commander 

US  Army  Missile  Command 

ATTN:  AMSMI-DA 

Redstone  Arsenal,  AL  35809 

1  Commander 

US  Army  Troop  Support  Command 
ATTN:  Sys  Anal  Ofc 
4300  Goodfellow  Blvd 
St.  Louis,  MO  63120 

1  Commander 

Rock  Island  Arsenal 
ATTN:  SARRI-LR-S 
Rock  Island,  IL  61201 

1  Commander 

US  Army  Tank-Automotive  Command 
ATTN:  AMSTA-RPL  (Tech  Lib) 
Warren,  MI  48090 

1  Commander 

US  Army  Armament  Command 
ATTN:  Sys  Anal  Ofc 
Rock  Island,  IL  61201 

1  Manager 

Defense  Logistics  Studies 
Information  Exchange 
US  Army  Logistics  Management 
Center 

ATTN:  Mr.  Leon  T.  Scarbrough 
Fort  Lee,  VA  23801 

Aberdeen  Proving  Ground 

Cdr,  USATECOM,  ATTN:  AMSTE-TS 
Ch,  Tech  Lib,  Bldg  305 
Dir,  BRL,  Bldg  328 


74 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PACE  (Whtn  Dmtm  Bnltnd) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT  NUMBER  2.  GOVT  ACCESSION  NO. 

TR-85 

3.  RECIPIENT'S  CATALOG  NUMBER 

4  TITLE  (and  Subtitle) 

A  Recursive  Estimation  Algorithm  for 
Identification  of  Dynamic  Systems 

5.  TYPE  OF  REPORT  &  PERIOD  COVERED 

6.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHORf*; 

Toney  R.  Perkins 

James  F.  Leathrum 

0.  CONTRACT  OR  GRANT  HUMBERft) 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADORESS 

Director 

U.  S.  Army  Materiel  Systems  Analysis  Agency 
Aberdeen  Proving  Ground,  Maryland  21005 

10.  PROGRAM  ELEMENT,  PROJECT,  TASK 
AREA  &  WORK  UNIT  NUMBERS 

RDTE  Project  No. 

1T765706M541 

It.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

Commander 

US  Army  Materiel  Command,  5001  Eisenhower  Avenue 
Alexandria.  VA 

12.  REPORT  DATE 

April  1974 

13.  NUMBER  OF  PAGES 

74 

14  MONITORING  AGENCY  NAME  &  ADDRESSf/f  diffarant  from  Controlling  Office) 

15  SECURITY  CLASS,  (of  thie  report) 

Unclassified 

15«.  DECLASS1  FI  CATION/ DOWNGRADING 

SCHEDULE 

16  DISTRIBUTION  STATEMENT  (of  thle  Report) 

Approved  for  public  release;  distribution  unlimited. 

17.  DISTRIBUTION  STATEMENT  (of  the  ebetrect  entered  In  Block  20.  If  different  from  Report ) 

Approved  for  public  release;  distribution  unlimited. 

18  SUPPLEMENTARY  NOTES 

19.  KEY  WORDS  (Continue  on  reverse  aide  If  n eceaaary  and  Identify  by  block  number) 

Parameter  identification,  recursive  estimation,  minimum  variance,  covariance, 
statistical  model,  prediction. 

20.  ABSTRACT  (Continue  on  reverse  aide  If  neceaeary  and  Identity  by  block  number) 

This  report  deals  with  the  identification  of  unknown  parameters  in  dynamic 
systems.  In  modeling  a  physical  system,  the  problem  of  identifying  the  dynamics 
of  a  system  are  often  encountered.  The  developed  algorithm  provides  a  tool  to 
model  all  or  parts  of  a  dynamic  system  using  input-output  data  sets  from  a  real 
system.  The  methodology  and  techniques  of  this  algorithm  are  based  upon  linear 
recursive  estimation  theory.  The  theoretical  foundation  and  the  pragmatics  of 
utilizing  the  ensemble  data  to  estimate  the  unknown  parameters  are  discussed  at 
length  in  the  development  of  the  algorithm.  (cont'd) 

DD  ,  j^M7J  1473  EDITION  OF  I  NOV  65  IS  OBSOLETE  UNCLASSIFIED 


A 


SECURITY  CL  ASS  I F I C  AT  ION  OF  THIS  PAGE  (When  Data  Entarad) 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  or  THIS  PAOEfPfc— 

Item  #20. 

As  an  application  of  the  algorithm,  experimental  data  from  a  man-in-the-loop 
simulation  is  used  to  estimate  the  parameters  of  a  single  axis  model  of  the 
gunner.  The  tracking  response  of  the  gunner  model  compare  favorably  with  data 
obtained  from  the  simulation.  The  differences  in  tracking  responses  are 
attributed  to  not  including  human  randomness  in  the  model. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE(T17>*n  Dmtm  Bntmr md) 


